Finite Elements Lab 2 Worksheet

Georges Limbert
nCATS, Faculty of Engineering and the Environment, University of Southampton

In [1]:
from IPython.core.display import HTML
css_file = '../ipython_notebook_styles/ngcmstyle.css'
HTML(open(css_file, "r").read())


Out[1]:

Aims

The aim of this computer laboratory is to gain a practical insight into the fundamental steps required to program a simple finite element code to solve a problem of elasticity. In doing so, you will learn:

  • How to recast the PDE problem into a variational form
  • To derive the formulation of a 2D finite element
  • To create a simple structured mesh generator to mesh a 2D domain
  • To create an assembly procedure for the stiffness matrix of a 2D structure
  • How to numerically solve the variational problem for a practical example

The numerical environment chosen for this laboratory is IPython.

Basic problem: linear elasticity

Let's consider a 2D solid body $\Omega $ at rest, subjected to body forces ${\bf{b}}$, prescribed displacements ${\bf{\bar d}}$ and tractions ${\bf{\bar t}}$. The problem of elasticity is defined as follows:

$$ \begin{align} {\text{div}}\sigma + {\bf{b}} & = {\bf{0}} & {\bf x}(x,y) &\text{ in } \Omega \\ {\bf{d}}(u,v) &= {\bf{\bar d}}(u,v) & {\bf x}(x,y) &\text{ on } \partial {\Omega _u} \\ {\bf{t}}({\bf x}) &= {\bf{\bar t}}({\bf x}) & {\bf x}(x,y) &\text{ on }\partial {\Omega_\sigma} \end{align} $$

$\sigma $ is the symmetric Cauchy stress tensor.

Question 1

Apply integration by parts to the PDE ${\text{div}}\sigma + {\bf{b}} = {\bf{0}}$ with the divergence theorem (relation between the flow (flux) of a vector field though a surface and the behavior inside the surface) recast it into its variational form. You will also take into account the conditions about the test function $\upsilon$ on the boundary where ${\bf{d}}$ is prescribed.

Constituitive relationship

To relate stress and strain one must define what is called a constitutive relationship. In other branches of physics other types of relationships would be used (e.g. relation between strain and temperature).

In the case of isotropic linear elasticity (Hooke's elasticity) the Cauchy stress tensor is related to the strain tensor $\varepsilon $ through the following relationship:

$$ \begin{equation} \sigma = \mathbb{C} : {\bf \varepsilon} \end{equation} $$

where $\mathbb{C}$ is the elasticity tensor.

Equation (2) is a tensor equation (valid in any coordinate system). If one decides to use a Cartesian coordinate system (2) can be written in index notation as:

$$ \begin{equation} {\left( {\bf \sigma} \right)_{ij}} = {\sigma _{ij}} = \left( \mathbb{C} : {\bf \varepsilon} \right)_{ij} = \mathbb{C}_{ijkl} \varepsilon_{kl} \quad \left\{ {{\text{2D}}:i,j,k,l:1 \dots 2; \quad {\text{3D}}:i,j,k,l:1 \dots 3} \right\} \end{equation} $$

Please note that we use the following convention: bold characters represent vectors, matrices or tensors.

Because the Cauchy, strain and elasticity tensors are symmetric and the material isotropic one can simplify the tensorial constitutive equation into a matrix equation involving vectors and asecond-order matrix.

$$ \begin{equation} \left\{ {\bf \sigma} \right\} = \left\{ \begin{array}{l} {\sigma _{xx}}\\ {\sigma _{yy}}\\ {\sigma _{xy}} \end{array} \right\} = \underbrace {\frac{E}{{1 - {\nu ^2}}}\left[ {\begin{array}{*{20}{c}} 1&\nu &0\\ \nu &1&0\\ 0&0&{\frac{{1 - \nu }}{2}} \end{array}} \right]}_{\left[ \mathbb{C} \right]}\underbrace {\left\{ \begin{array}{l} {\varepsilon _{xx}}\\ {\varepsilon _{yy}}\\ 2{\varepsilon _{xy}} \end{array} \right\}}_{\left\{ {\bf \varepsilon} \right\}} \end{equation} $$

$\left\{ {E,\nu } \right\}$ are material properties, respectively the Young's modulus and Poisson's ratio which characterise the stiffness and compressibility of the material.

Let’s consider a 4-noded bi-linear interpolation element (Figure 1).

The shape functions chosen here are simple linear Lagrange polynomials but other forms could be used (quadratic Lagrange, Hermite, NURBS, etc):

\begin{equation} {N_i} = \frac{1}{4}(1 + {\xi_i} \xi )(1 + {\eta _i}\eta ) \quad {\xi _i},{\eta _i} = - 1,1 \end{equation}

There are as many shape functions as nodes in the element, i.e. 4:

\begin{equation} \left\{ \begin{array}{rl} N_1 &= \frac{1}{4}\left( {1 - \xi } \right)\left( {1 - \eta } \right)\\ N_2 &= \frac{1}{4}\left( {1 + \xi } \right)\left( {1 - \eta } \right)\\ N_3 &= \frac{1}{4}\left( {1 + \xi } \right)\left( {1 + \eta } \right)\\ N_4 &= \frac{1}{4}\left( {1 - \xi } \right)\left( {1 + \eta } \right) \end{array} \right\} \end{equation}

We are using an isoparametric formulation which means that the position vector and displacement vector are interpolated using the same order of shape (or interpolation) functions.

Interpolation of the position vector (discretisation)

$$ \begin{equation} {\bf{x}}\left( {{\xi ^e},t} \right) = {{\bf{x}}_1}(t){N_1}\left( {{\xi ^e}} \right) + {{\bf{x}}_2}(t){N_2}\left( {{\xi ^e}} \right) + {{\bf{x}}_3}(t){N_3}\left( {{\xi ^e}} \right) + {{\bf{x}}_4}(t){N_4}\left( {{\xi ^e}} \right) \end{equation} $$

where

$$ \begin{equation} {{\bf{x}}_A}(t) = {\rm{position}}\,(node\,A) = \left\{ \begin{array}{l} {x_A}\\ {y_A} \end{array} \right\} \end{equation} $$$$ \begin{equation} {\bf{x}}\left( {{\xi ^e},t} \right) = {\left\{ {\begin{array}{*{20}{c}} {x(t)}\\ {y(t)} \end{array}} \right\}_{node\,1}}{N_1}\left( {{\xi ^e}} \right) + {\left\{ {\begin{array}{*{20}{c}} {x(t)}\\ {y(t)} \end{array}} \right\}_{node\,2}}{N_2}\left( {{\xi ^e}} \right) + {\left\{ {\begin{array}{*{20}{c}} {x(t)}\\ {y(t)} \end{array}} \right\}_{node\,3}}{N_3}\left( {{\xi ^e}} \right) + {\left\{ {\begin{array}{*{20}{c}} {x(t)}\\ {y(t)} \end{array}} \right\}_{node\,4}}{N_4}\left( {{\xi ^e}} \right) \end{equation} $$

Interpolation of the displacement vector (discretisation)

$$ \begin{equation} {\bf{d}}\left( {{\xi ^e},t} \right) = {{\bf{d}}_1}(t){N_1}\left( {{\xi ^e}} \right) + {{\bf{d}}_2}(t){N_2}\left( {{\xi ^e}} \right) + {{\bf{d}}_3}(t){N_3}\left( {{\xi ^e}} \right) + {{\bf{d}}_4}(t){N_4}\left( {{\xi ^e}} \right) \end{equation} $$

where

$$ \begin{equation} {{\bf{d}}_A}(t) = {\rm{displacement}}\,(node\,A) = \left\{ \begin{array}{l} {u_A}\\ {v_A} \end{array} \right\} \end{equation} $$

Degrees of freedom

We can arrange all the degrees of freedom of the element in one vector:

$$ \begin{equation} {\bf{d}} = {\{ {u_1},{v_1},{u_2},{v_2},{u_3},{v_3},{u_4},{v_4}\} ^T} \end{equation} $$

The discretised versions of the position and displacement vector can be rewritten in matrix form:

$$ \begin{equation} {\bf{x}}\left( {{\xi ^e},t} \right) = {N_I}\left( {{\xi ^e}} \right){{\bf{x}}_I}(t) = {{\bf{N}}^T}\left( {{\xi ^e}} \right){{\bf{x}}_I}(t) \end{equation} $$$$ \begin{equation} {\bf{d}}\left( {{\xi ^e},t} \right) = {N_I}\left( {{\xi ^e}} \right){{\bf{d}}_I}(t) = {{\bf{N}}^T}\left( {{\xi ^e}} \right){{\bf{d}}_I}(t) \end{equation} $$

where

$$ \begin{equation} {{\bf{N}}^T}\left( {{\xi ^e}} \right) = \left\{ {{N_1}({\xi ^e}),{N_2}({\xi ^e}),{N_3}({\xi ^e}),{N_4}({\xi ^e})} \right\} \end{equation} $$

The test function used for the formulation of the weak form is also discretised:

$$ \begin{equation} \upsilon \left( {{\xi ^e},t} \right) = {N_I}\left( {{\xi ^e}} \right){\upsilon _I}(t) = {{\bf{N}}^T}\left( {{\xi ^e}} \right){\upsilon _I}(t) \end{equation} $$

Question 2

Inject the discretised version of the displacement into the weak form established in Question 1. You will have to express the strain matrix using the ${\bf B}$ matrix (derivatives of shape functions). You should obtain 3 integral matrix/vector forms: the element stiffness matrix, the nodal load vector due to body forces and the nodal load vector due to boundary forces.

Integration

The integral expressions obtained in Question 2 will have to be numerically integrated using, in our case, Gauss integration. We have integrals defined over the element domain ${\Omega ^e}$ and its boundary $\partial {\Omega ^e}$.

There are well known integration formulas over parametric domains that are exact for polynomials like the Lagrange polynomials used in the shape functions. The integral of a function is transformed into the sum of weighted values that the function takes at specific points, the so-called integration points or Gauss points.

$$ \begin{equation} \int\limits_ {f(\xi )\,\,} d \xi = {w_i}f({\xi _i}) \end{equation} $$

These integration formulas apply to the parent element (Figure 1) which is the element in the parametric domain $ = \{ \{ - 1,1\} ,\{ - 1,1\} \} $ so it is necessary to map back the parent element domain to the current (deformed) domain. This is a simple change of variables:

$$ \begin{equation} d{\Omega _e} = dx\,dy = \det ({\bf{J}})\,d\xi\, d\eta \end{equation} $$$$ \begin{equation} \int\limits_{{\Omega^{e}}} {f({\bf{x}})\,\,} d\Omega = \int\limits_ {f({\bf{x}}({\xi ^e}))\,\underbrace {\det \left( {\frac{{\partial {\bf{x}}}}{{\partial {\xi ^e}}}} \right)}_{{J_{{\xi ^e}}} = \det ({\bf{J}})}\,} d \xi = {w_i}{J_{{\xi ^e}}}({\xi ^e}_i)f({\xi ^e}_i) \end{equation} $$

Here, we will use a 2 points Gauss integration rule: i.e. 4 integration points for the quadrilateral element:

$$ \begin{align} {{\bf{G}}_1} &= \left\{ { - \frac{1}{{\sqrt 3 }}, - \frac{1}{{\sqrt 3 }}} \right\} & {w_1} &= 1\\ {{\bf{G}}_2} &= \left\{ {\frac{1}{{\sqrt 3 }}, - \frac{1}{{\sqrt 3 }}} \right\} & {w_2} &= 1\\ {{\bf{G}}_3} &= \left\{ {\frac{1}{{\sqrt 3 }},\frac{1}{{\sqrt 3 }}} \right\} & {w_3} &= 1\\ {{\bf{G}}_4} &= \left\{ { - \frac{1}{{\sqrt 3 }},\frac{1}{{\sqrt 3 }}} \right\} & {w_4} &= 1 \end{align} $$

Question 3

  1. Inject the integration formula into the integral matrix equations obtained in Question 2.
  2. Write a Python programme to:

    1. calculate the stiffness matrix of the element and nodal force vector;
    2. generate a structured grid over a 2D rectangle domain (mesh generator);
    3. specify a uniform force at each node of one side of the rectangle domain;
    4. specify a uniform displacement at each node of one side of the rectangle domain;
    5. assemble the global stiffness matrix and force vector of the whole rectangular structure.

Question 4

Using the Python programme from question 3, solve the following problem.

Consider a 2D rectangular solid domain (length $L_x = 1m$; height $L_y = 0.5m$) with material properties $\{ E = 210{\rm{ GPa}};\,\,\,\nu = 0.3\} $ and sides labelled (counter-clockwise, starting from the left vertical side) as A, B, C and D.

The boundary conditions are as follows:

\begin{equation} u({\text{side}}\,{\text{A}}) = v({\text{side}}\,{\text{B}}) = 0 \end{equation}

The load condition is:

\begin{equation} {\bf{F}}({\text{side C}})[Newton] = \left\{ \begin{array}{l} {F_x}({\text{side C}})\\ {F_y}({\text{side C}}) \end{array} \right\} = \left\{ \begin{array}{l} {10^6}\\ 0 \end{array} \right\} \end{equation}

The boundary and loading conditions are homogeneous and therefore do not introduce shear coupling. The Cauchy stress tensor has only one non-null component, ${\sigma _{xx}}$.

NB: the total load of 10$^{6}$ N must be uniformly distributed over all the nodes of side C.

  1. Solve the boundary-value problem of elasticity described above to calculate the displacement of faces D and C. You will consider 2 mesh sizes: 1 for which the whole structure is meshed with one element only and one size of your choice.
  2. Compare the solutions to the exact analytical solutions given by ${\sigma _{xx}} = {F_x}/{L_y}$.

The strains can be calculated as

$$ \begin{align} {\varepsilon _{xx}} = {\sigma _{xx}}/E = {F_x}/E{L_y} &= - {\varepsilon _{yy}}/\nu \\ {\varepsilon _{xy}} &= 0 \end{align} $$

and the displacements as

$$ \begin{align} u &= {\varepsilon _{xx}}{L_x} \\ v &= {\varepsilon _{yy}}{L_y}. \end{align} $$

Do the results surprise you? Why?